Thalamocortical Connectivity: Atlas Characteristics and Tract Anatomy
Thalamocortical Autotrack Atlas Anatomical Characteristics
Glasser region names/labels
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv") #mesulam assignments for glasser parcels
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)Spatial permutations of the glasser parcellation for spin testing
#https://github.com/frantisekvasa/rotate_parcellation/tree/master, cloned 8/5/2023
source("/cbica/projects/thalamocortical_development/software/rotate_parcellation/R/rotate.parcellation.R")
source("/cbica/projects/thalamocortical_development/software/rotate_parcellation/R/perm.sphere.p.R")
glasser.coords <- read.table("/cbica/projects/thalamocortical_development/software/rotate_parcellation/sphere_HCP.txt", header=F) #The spherical centroids of the multimodal HCP parcellation on the freesurfer sphere; the 360 regions are originally ordered here as left hemisphere -> right hemisphere. The glasser region list and the S-A axis are ordered right hemisphere -> left hemisphere however so let's update the ordering here for spin tests to ensure correspondence
glasser.coords <- rbind(glasser.coords[181:360,], glasser.coords[1:180,])
perm.id.full <- rotate.parcellation(coord.l = as.matrix(glasser.coords[181:360,]), coord.r = as.matrix(glasser.coords[1:180,]), nrot = 10000, method = "vasa") #rotate the glasser parcellation 10,000 times on the freesurfer sphere to generate spatial nulls for spin-based permutation significance testing
saveRDS(perm.id.full, "/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds")perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins for spin testing (spherical permutation based p-values). Parcel order is right hemisphere -> left hemisphere
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinningGlasser parcel anatomy: surface area and sulcal depth measures
glasserparcel.anatomy <- read.csv("/cbica/projects/thalamocortical_development/Templates/fsaverage/stats/glasserparcel_anatomy.csv") %>% select(StructName, Area_mm2, Mean) %>% filter(StructName != "???") %>% set_names("orig_parcelname", "surface_area", "sulcal_depth") #parcel name, total surface area, and mean sulcal depth for 360 glasser regions, calculated in /parcel_anatomy/glasser_parcel_anatomy bash and python scripts
glasserparcel.anatomy <- merge(glasserparcel.anatomy, glasser.assignments, by = "orig_parcelname", sort = F)Thalamocortical connection atlas tract list
tc.autotrack.tracts <- read.csv("/cbica/projects/thalamocortical_development/software/thalamocortical_autotrack_template/dsi-studio/atlas/ICBM152_adult/ICBM152_adult.tt.gz.txt", header = F) #238 connections
colnames(tc.autotrack.tracts) <- c("tract")
tc.autotrack.tracts$tract <- gsub("-", "_", tc.autotrack.tracts$tract) #no dasher -, no prancer, no vixenSpin-based enrichment testing function
#A function to test whether an anatomical measure of interest is significant larger or smaller in cortical regions without thalamic connections included in the atlas, based on a spin-based null distribution of the measure
source("/cbica/projects/thalamocortical_development/code/thalamocortical_development/results/tract_anatomy/tractometryatlas_enrichment.R")Cortical Surface Coverage
Percent of the cortical surface reached by thalamocortical atlas connections
#total area of the cortical surface (sum of the area of all glasser parcels)
total.surfacearea <- sum(glasserparcel.anatomy$surface_area)
#area of the cortical surface with reconstructed thalamus connections
tract.surfacearea <- sum(glasserparcel.anatomy[glasserparcel.anatomy$tract %in% tc.autotrack.tracts$tract,]$surface_area)
#percent of surface area with thalamus-cortex structural connections
round((tract.surfacearea/total.surfacearea)*100,2)## [1] 76.82
Anatomical Correlates of Reconstructed Connections
Surface area of parcels with versus without reconstructed thalamocortical connections
glasserparcel.anatomy$reconstructed <- ifelse(glasserparcel.anatomy$tract %in% tc.autotrack.tracts$tract, "Present", "Absent")
glasserparcel.anatomy$reconstructed <- factor(glasserparcel.anatomy$reconstructed, ordered = T, levels = c("Present", "Absent"))glasserparcel.anatomy <- glasserparcel.anatomy %>% mutate(surface_area_z = scale(surface_area, center = T, scale = T)) #z-score surface areaggplot(glasserparcel.anatomy, aes(x = reconstructed, y = surface_area_z)) +
geom_point(aes(color = reconstructed, fill = reconstructed), position = position_jitter(width = .4), size = 0.07) +
geom_boxplot(aes(color = reconstructed), fill = alpha("white", 0.75), size = 0.15, outlier.shape = NA) +
theme_minimal() +
scale_color_manual(values = c("#48457F", "#FCAB6A")) +
labs(x="\n", y="surface area\n") +
theme(legend.position = "none") +
theme(
panel.grid.major = element_line(linewidth = 0.2),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1/Atlasconnections_surfacearea.pdf", device = "pdf", dpi = 500, width = 1.3 , height = 1.2)## [1] 0
Sulcal depth of parcels with versus without reconstructed thalamocortical connections
#calculate percentage of maximal sulcal depth
#in freesurfer, gyri have negative values and sulci have positive values as they reflect distance from a mid-surface of 0. Let's turn this into percent of maximum sulcal depth
total_sulcalrange <- max(glasserparcel.anatomy$sulcal_depth) - min(glasserparcel.anatomy$sulcal_depth) #span of sulcal depth data
distance_from_gyralcrown <- glasserparcel.anatomy$sulcal_depth - (min(glasserparcel.anatomy$sulcal_depth)) #relative distance from gyral crown
glasserparcel.anatomy$percent_sulcaldepth <- (distance_from_gyralcrown/total_sulcalrange)*100 #percent of total depthggplot(glasserparcel.anatomy, aes(x = reconstructed, y = percent_sulcaldepth)) +
geom_point(aes(color = reconstructed, fill = reconstructed), position = position_jitter(width = .4), size = 0.07) +
geom_boxplot(aes(color = reconstructed), fill = alpha("white", 0.75), size = 0.25, outlier.shape = NA) +
theme_minimal() +
scale_color_manual(values = c("#48457F", "#FCAB6A")) +
labs(x="\n", y="Percentage sulcal depth\n") +
theme(legend.position = "none") +
theme(
panel.grid.major = element_line(linewidth = 0.2),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1/Atlasconnections_sulcaldepth.pdf", device = "pdf", dpi = 500, width = 1.3 , height = 1.2)## [1] 0.0348
Thalamic Nuclei-specific Cortical Connection Zones
#Number of voxels occupied by each thalamocortical connection in each thalamic nucleus, generated by /thalamocortical_structuralconnectivity/template/thalamocortical_connection_nuclei scripts and based on the Morel-guided MR atlas from Saranathan et al., 2021, Scientific Data
nucleus.connection.profiles <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/template/HCP-MMP_thalamicnuclei_terminationzones.csv") %>% select(tract, label2, label4, label5, label6, label7, label8, label9, label10, label11, label12) #labels 1 and 3 are not used *"unlabeled") in this atlas
colnames(nucleus.connection.profiles) <- c("tract", "AV", "VA", "VLa", "VLp", "VPL", "Pul", "LGN", "MGN", "CM", "MD") #label-nucleus assignments from /cbica/projects/thalamocortical_development/Maps/thalamicnuclei_atlas_saranathan/CustomAtlas.ctbl
connection.list <- gsub("-", "_", nucleus.connection.profiles$tract)#A function to plot cortical regions that show strongest connectivity to selected thalamic nuclei
nucleus.connection.plot <- function(nuclei.names, n.regions, cortex.color){
nucleus.connection.profiles %>% select(all_of(nuclei.names)) %>% #select thalamic nuclei
rowSums %>% as.data.frame() %>% set_names("connectivity.strength") %>% #sum connectivity counts
mutate(tract = connection.list) %>% merge(., glasser.regions, by = "tract") %>%
slice_max(order_by = connectivity.strength, n = n.regions, with_ties = F) %>% #select top connected cortical regions
mutate(top.connected = "yes") %>%
select(label, top.connected) %>% mutate(cortex = "cortex") %>%
rbind(., data.frame(label = "rh_???", top.connected = "no", cortex = "medialwall")) %>%
rbind(., data.frame(label = "lh_???", top.connected = "no", cortex = "medialwall")) %>%
filter(grepl("lh_", label)) %>%
ggplot(.) +
geom_brain(atlas = glasser, hemi = "left", mapping = aes(fill = top.connected, colour=I("black"), size=I(.05))) +
theme_void() +
scale_fill_manual(values = c("white", cortex.color), na.value = "white") +
new_scale_fill() +
geom_brain(atlas = glasser, hemi = "left", mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c(alpha("white", 0), "grey"), na.value = "white") +
theme(legend.position = "none")
}VLa and VLp
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/VLaVLp_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
VPL
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/VPL_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
Pulvinar, LGN
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/PulLGN_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
CM
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/CM_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
MD
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/MD_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
AV
## merging atlas and data by 'label'
## merging atlas and data by 'label'
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure1_supplement/AV_connectivity.pdf", device = "pdf", dpi = 500, width = 2.5, height = 3.5)## merging atlas and data by 'label'
## merging atlas and data by 'label'
Thalamocortical Connectivity Cortical and Thalamic Gradients
S-A axis ranks
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by="orig_parcelname", sort = F)Dataset-specific final tract lists
#Tract lists generated by /tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv")
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv")Function to calculate tract-level means per dataset
tract_average_measures <- function(measure, atlas, dataset){
input.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset)) %>% select(contains("autotrack")) #nrow subjects by ncol tracts for mean calculation
tract.means <- colMeans(input.df, na.rm = TRUE) %>% as.data.frame %>% set_names(sprintf("mean_%s", measure)) %>% mutate(tract = row.names(.))
tract.means <- merge(tract.means, S.A.axis, by = "tract")
return(tract.means)
}Thalamocortical Connection Fractional Anisotropy
FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #FA for included tracts and participants
FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") Calculate mean FA for each connection in PNC and HCPD datasets
tract.meanFA.pnc <- tract_average_measures(measure = "FA", atlas = "glasser", dataset = "pnc")
tract.meanFA.pnc <- tract.meanFA.pnc[tract.meanFA.pnc$tract %in% tractlist.pnc$tract,] #connection-level mean FA for final PNC tract list
write.csv(tract.meanFA.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_FA_primary.csv", row.names = F, quote = F)
tract.meanFA.hcpd <- tract_average_measures(measure = "FA", atlas = "glasser", dataset = "hcpd")
tract.meanFA.hcpd <- tract.meanFA.hcpd[tract.meanFA.hcpd$tract %in% tractlist.hcpd$tract,] #connection-level mean FA for final HCPD tract list
write.csv(tract.meanFA.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_FA_primary.csv", row.names = F, quote = F)Correlation of connection-level average FA between datasets
Pearson’s r
tract.meanFA.merged <- merge((tract.meanFA.pnc %>% select(tract, mean_FA)), (tract.meanFA.hcpd %>% select(tract, mean_FA)), by = "tract") %>% set_names("tract", "mean_FA.pnc", "mean_FA.hcpd")
cor.test(tract.meanFA.merged$mean_FA.pnc, tract.meanFA.merged$mean_FA.hcpd)##
## Pearson's product-moment correlation
##
## data: tract.meanFA.merged$mean_FA.pnc and tract.meanFA.merged$mean_FA.hcpd
## t = 61.101, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9635117 0.9783814
## sample estimates:
## cor
## 0.9719002
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, tract.meanFA.merged, by = "tract") #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$mean_FA.pnc, y = spin.data$mean_FA.hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0
Correlation plot
ggplot(tract.meanFA.merged, aes(x = mean_FA.pnc, y = mean_FA.hcpd)) +
geom_point(color = c("#571766"), size = 0.5) +
geom_smooth(method = "lm", color = c("goldenrod1"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(0.2, 0.525)) +
scale_y_continuous(limits = c(0.2, 0.525)) +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/FA_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.4)Thalamocortical connection FA organization along the S-A axis
Spearman’s correlations
PNC
##
## Spearman's rank correlation rho
##
## data: tract.meanFA.pnc$mean_FA and tract.meanFA.pnc$SA.axis
## S = 2633482, p-value = 4.941e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4248717
HCPD
##
## Spearman's rank correlation rho
##
## data: tract.meanFA.hcpd$mean_FA and tract.meanFA.hcpd$SA.axis
## S = 2877906, p-value = 4.757e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4192293
Spatial permutation (spin) based p-value
PNC
spin.data <- left_join(spin.parcels, tract.meanFA.pnc, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_FA, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.045
HCPD
spin.data <- left_join(spin.parcels, tract.meanFA.hcpd, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_FA, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.06135
Thalamocortical FA by S-A axis plot
PNC
ggplot(tract.meanFA.pnc, aes(x = SA.axis, y = mean_FA)) +
geom_point(aes(color = SA.axis), size = 0.5) +
geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) +
theme_classic() +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
labs(x="\nSensorimotor-association axis rank", y="Connection FA\n") +
theme(legend.position = "none") +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/FA_SArank_PNC.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)HCPD
ggplot(tract.meanFA.hcpd, aes(x = SA.axis, y = mean_FA)) +
geom_point(aes(color = SA.axis), size = 0.5) +
geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) +
theme_classic() +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
labs(x="\nSensorimotor-association axis rank", y="Connection FA\n") +
theme(legend.position = "none") +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical Connection Core-Matrix Gradient (CPt) Values
Read in connection-specific core-matrix gradient values (based on the connection’s thalamic endpoints and the thalamic calbindin-parvalbumin CPt gradient) and get data for just final study sample and included tracts
extract_measure <- function(metric.df, participants.df, measure){
measure.df <- metric.df %>% select(rbcid, tract, all_of(measure)) #extract measure of interest from long df
measure.df <- measure.df %>% arrange(tract) #tracts in alphabetical order
measure.df <- measure.df %>% pivot_wider(names_from = "tract", values_from = all_of(measure)) #long to wide format, all participants and tracts
measure.df <- measure.df %>% arrange(rbcid) #rbcids in numerical order
colnames(measure.df) <- gsub("-", "_", colnames(measure.df))
measure.df <- measure.df[measure.df$rbcid %in% participants.df$rbcid,] #only final sample of participants now
return(measure.df)
}PNC
#Create connection CPt value spreadsheet for final sample of PNC participants
CPt.glasser.pnc.long <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_thalamicconnection_CPt.csv") #thalamic connection CPt gradient values, generated by /tract_measures/CPtgradient_values/lthalamocortical_calculate_CPtvalues.sh and /tract_measures/CPtgradient_values/thalamocortical_CPtvalues.py
participants.pnc <- read.csv("/cbica/projects/thalamocortical_development/sample_info/PNC/PNC_thalamocortical_finalsample_N1145.csv") #PNC final participant list
CPt.glasser.pnc <- extract_measure(metric.df = CPt.glasser.pnc.long, participants.df = participants.pnc, measure = "CPt")#Don't analyze CPt values for connections with <5 streamlines
streamlinecount.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_streamlinecount_finalsample.csv")
identical(streamlinecount.glasser.pnc$rbcid, CPt.glasser.pnc$rbcid) ## [1] TRUE
streamlinecount.glasser.pnc <- streamlinecount.glasser.pnc %>% select(contains("autotrack")) #1145 subjects by 238 tracts
CPt.glasser.pnc <- CPt.glasser.pnc %>% select(contains("autotrack")) #1145 subjects by 238 tracts
identical(names(streamlinecount.glasser.pnc), names(CPt.glasser.pnc))## [1] TRUE
HCPD
#Create connection CPt value spreadsheet for final sample of HCPD participants
CPt.glasser.hcpd.long <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_thalamicconnection_CPt.csv")
participants.hcpd <- read.csv("/cbica/projects/thalamocortical_development/sample_info/HCPD/HCPD_thalamocortical_finalsample_N572.csv") #HCPD final participant list
CPt.glasser.hcpd <- extract_measure(metric.df = CPt.glasser.hcpd.long, participants.df = participants.hcpd, measure = "CPt")#Don't analyze CPt values for connections with <5 streamlines
streamlinecount.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_streamlinecount_finalsample.csv")
identical(streamlinecount.glasser.hcpd$rbcid, CPt.glasser.hcpd$rbcid)## [1] TRUE
streamlinecount.glasser.hcpd <- streamlinecount.glasser.hcpd %>% select(contains("autotrack")) #572 subjects by 238 tracts
CPt.glasser.hcpd <- CPt.glasser.hcpd %>% select(contains("autotrack")) #572 subjects by 238 tracts
identical(names(streamlinecount.glasser.hcpd), names(CPt.glasser.hcpd)) ## [1] TRUE
Calculate thalamic connection area CPt value for each tract in PNC and HCPD datasets
tract.meanCPt.pnc <- tract_average_measures(measure = "CPt", atlas = "glasser", dataset = "pnc")
tract.meanCPt.pnc <- tract.meanCPt.pnc[tract.meanCPt.pnc$tract %in% tractlist.pnc$tract,]
write.csv(tract.meanCPt.pnc, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractaverage_CPt_primary.csv", row.names = F, quote = F)
tract.meanCPt.hcpd <- tract_average_measures(measure = "CPt", atlas = "glasser", dataset = "hcpd")
tract.meanCPt.hcpd <- tract.meanCPt.hcpd[tract.meanCPt.hcpd$tract %in% tractlist.hcpd$tract,]
write.csv(tract.meanCPt.hcpd, "/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractaverage_CPt_primary.csv", row.names = F, quote = F)Correlation of connection-level average CPt between datasets
Pearson’s r
tract.meanCPt.merged <- merge((tract.meanCPt.pnc %>% select(tract, mean_CPt)), (tract.meanCPt.hcpd %>% select(tract, mean_CPt)), by = "tract") %>% set_names("tract", "mean_CPt.pnc", "mean_CPt.hcpd")
cor.test(tract.meanCPt.merged$mean_CPt.pnc, tract.meanCPt.merged$mean_CPt.hcpd)##
## Pearson's product-moment correlation
##
## data: tract.meanCPt.merged$mean_CPt.pnc and tract.meanCPt.merged$mean_CPt.hcpd
## t = 93.485, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9839914 0.9905554
## sample estimates:
## cor
## 0.9877012
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, tract.meanCPt.merged, by = "tract")
perm.sphere.p(x = spin.data$mean_CPt.pnc, y = spin.data$mean_CPt.hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0
Correlation plot
ggplot(tract.meanCPt.merged, aes(x = mean_CPt.pnc, y = mean_CPt.hcpd)) +
geom_point(color = c("#571766"), size = 0.5) +
geom_smooth(method = "lm", color = c("goldenrod1"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/CPt_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.4)Thalamocortical connection CPt organization along the S-A axis
Spearman’s correlations
PNC
##
## Spearman's rank correlation rho
##
## data: tract.meanCPt.pnc$mean_CPt and tract.meanCPt.pnc$SA.axis
## S = 781550, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5771346
HCPD
##
## Spearman's rank correlation rho
##
## data: tract.meanCPt.hcpd$mean_CPt and tract.meanCPt.hcpd$SA.axis
## S = 1015854, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4990352
Spatial permutation (spin) based p-value
PNC
spin.data <- left_join(spin.parcels, tract.meanCPt.pnc, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_CPt, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.00835
HCPD
spin.data <- left_join(spin.parcels, tract.meanCPt.hcpd, by = c("orig_parcelname", "label", "tract"))
perm.sphere.p(x = spin.data$mean_CPt, y = spin.data$SA.axis, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.04045
Thalamocortical FA by S-A axis plot
PNC
ggplot(tract.meanCPt.pnc, aes(x = SA.axis, y = mean_CPt)) +
geom_point(aes(color = SA.axis), size = 0.5) +
geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) +
theme_classic() +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
labs(x="\nSensorimotor-association axis rank", y="Connection CPt\n") +
theme(legend.position = "none") +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure3/CPt_SArank_PNC.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.4)HCPD
ggplot(tract.meanCPt.hcpd, aes(x = SA.axis, y = mean_CPt)) +
geom_point(aes(color = SA.axis), size = 0.5) +
geom_smooth(method = "lm", linewidth = 0.3, color = c("gray25")) +
theme_classic() +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, limits = c(1, 360), midpoint = 180) +
labs(x="\nSensorimotor-association axis rank", y="Connection CPt\n") +
theme(legend.position = "none") +
theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))